Time Series Analysis of U.S. Monthly Tourist Arrivals and Spending

I. Introduction

The economic impact of international tourism in the United States is substantial, with global travelers contributing billions annually to the U.S. economy. Before the COVID-19 pandemic, international visitors spent approximately $233.5 billion in 2019 alone, which averages to nearly $640 million daily (International Trade Administration, n.d.). With U.S. tourism accounting for 14.5% of total global spending on international travel (Zaheer, 2023), understanding the trends and drivers behind these expenditures is vital for sustaining this economic sector.

While tourist spending is influenced by various socioeconomic factors, one of the most reliable indicators of such expenditures is the volume of international arrivals, commonly referred to as “Tourist Arrivals.” Fluctuations in arrivals often correlate closely with factors like currency exchange rates, seasonal trends, and global geopolitical stability (Hinton, 2024). In an attempt to explore the dynamics between tourist arrivals and spending, this project seeks to find the best predictive models of both metrics. Leveraging upon the function of time and SARIMA models, this study aims to provide a robust and comprehensive understanding of U.S. tourism’s economic contributions, as well as a tool for predicting future trends in international arrivals and spending.

II. Methods

1. Working Data

a. Data Acquisition

Arrivals <- read.csv("Monthly_Tourist_Arrivals_R.csv")
Exports <- read.csv("Monthly_Exports_R.csv")

Arrivals$Total <- as.numeric(Arrivals$Total)
Arrivals$Time <- floor_date(as.Date(Arrivals$Time), "month")
Exports$Time <- as.Date(Exports$Time)

Masterdata0 <- merge(Arrivals, Exports, by = "Time", all = TRUE)

For this project, we sourced two datasets from the International Trade Administration (ITA) website, both available in Excel format. The first dataset, “Monthly Exports Imports Balance,” was retrieved from the ITA’s International Travel Receipts and Payments Program. It provides monthly tourism-related trade data from January 1999 to July 2024, including metrics like total travel exports (receipts, or tourist spending) and imports (payments, or US citizen spending abroad), which are further subdivided into categories such as Travel Spending, Medical/ Education/ Workers Spending, and Passenger Fare Receipts. All figures are in millions of dollars and are seasonally adjusted.

The second dataset, “Monthly Arrivals 2000 to Present – Country of Residence (COR),” was acquired from the ITA’s I-94 Arrivals Program. This file contains monthly data on the number of tourists arriving in the U.S. from January 2000 to July 2024, with a breakdown by world region (e.g., Central America, Western Europe) and individual countries of residence. This dataset also includes year-over-year changes for each category. No specific filters or exclusions were applied to the data before downloading.

b. Potential Biases

The datasets used in this study inherently exclude certain categories of visitors. For instance, the International Visitor Arrivals Program - ADIS I-94 does not count travelers who do not stay overnight in the U.S. or those who visit on certain visa types. This exclusion may introduce some problems, as short-term visitors, particularly those near borders or in airports, may still contribute significantly to tourism-related spending but are not included in official tourist arrival figures. Another potential bias is related to the seasonality of tourist arrivals, which can be influenced by holidays and business cycles. Significant anomalies like the COVID-19 pandemic have also impacted arrival and spending patterns.

c. Data Wrangling

To facilitate our analysis, we simplified both datasets by focusing on aggregated monthly data. For the “Monthly Exports Imports Balance” file, we retained only the total tourism-related exports (tourist spending). Similarly, for the “Monthly Arrivals” file, we retained only the total arrivals and subtotal information by world region, removing individual country data to streamline our analysis. The figure below shows a quick comparison of the two main time series remaining.

plot1 <- ggplot(Masterdata0, aes(x = Time, y = Total)) +
  geom_line(color = "darkblue", size = 1) +
  ggtitle("US Monthly Tourist Arrivals (2000-2024)") +
  xlab("Month") +
  ylab("Tourist Arrivals") +
  scale_x_date(date_labels = "%b %Y", date_breaks = "36 months")+
  theme_light()+
  scale_y_continuous(labels = scales::label_number(scale_cut = cut_short_scale()))

plot2 <- ggplot(Masterdata0, aes(x = Time, y = Exports)) +
  geom_line(color = "darkred", size = 1) +
  ggtitle("US Monthly Tourist Exports, seasonally adjusted (2000-2024)") +
  xlab("Month") +
  ylab("Export Values (million $)") +
  scale_x_date(date_labels = "%b %Y", date_breaks = "36 months")+
  theme_light()+
  scale_y_continuous(labels = scales::label_number(scale_cut = cut_short_scale()))

grid.arrange(plot1, plot2, nrow=2)

Given the unique impact of COVID-19 on travel and spending, we intentionally excluded data post-January 2020. This choice allows us to analyze trends under more typical conditions, thus providing a more stable basis for examining seasonal and trend-based patterns in tourist arrivals and spending from 2000 to 2019. The figure below illustrates such exclusion of data.

Masterdata <-Masterdata0 %>%
  filter(Time < as.Date("2020-01-01"))%>%
  select("Time","Total","Exports") %>%
  mutate("Month" = month(Time)) %>%
  mutate("Year" = year(Time))

Arrivals_ts = ts(Masterdata$Total,start= c(2000,1), frequency = 12)
Exports_ts = ts(Masterdata$Exports,start= c(2000,1), frequency = 12)

tstools::tsplot(list("Arrivals (millions), left axis"=Arrivals_ts/1000000),
                tsr = list("Exports (billions $), right axis" =Exports_ts/1000), 
                theme = init_tsplot_theme(line_colors = c("darkblue", "darkred")),
                plot_title = "US Monthly Tourist Arrivals & Spending (2000-2019)")

To better facilitate the cross-correlation part between the two time series, we de-seasonalize the Tourist Arrivals time series by Seasonal Index Approach (more information in the Apendix). This is because we believe seasonally adjusted values should work better together, particularly when compared to one seasonlly adjusted and one normal. The seasonally adjusted values are plotted below.

# STL Decomposition
stl_arrivals <- stl(Arrivals_ts, s.window = 12)
Arrivals_sa_stl <- seasadj(stl_arrivals)

tstools::tsplot(list("Arrivals (millions), left axis"=Arrivals_sa_stl/1000000),
                tsr = list("Exports (billions $), right axis" =Exports_ts/1000), 
                theme = init_tsplot_theme(line_colors = c("darkblue", "darkred")),
                plot_title = "US Monthly Tourist Arrivals & Spending, seasonally adjusted (STL) (2000-2019)")

# SAAR Adjustment
Arrivals_SA <- Masterdata %>%
  select(c(Time, Total, Year, Month)) %>%
  group_by(Year) %>%
  mutate("Annual_Average" = mean(Total))%>%
  mutate("Yearly_Percentage" = Total/Annual_Average)

Arrivals_SA <- Arrivals_SA %>%
  group_by(Month)  %>%
  mutate("Seasonal_Index" = mean(Yearly_Percentage)) %>%
  mutate("SA_Arrivals" = Total/Seasonal_Index)
  
Arrivals_sa_saar <- ts(Arrivals_SA$SA_Arrivals,start= c(2000,1), frequency = 12)


tstools::tsplot(list("Arrivals (millions), left axis"=Arrivals_sa_saar/1000000),
                tsr = list("Exports (billions $), right axis" =Exports_ts/1000), 
                theme = init_tsplot_theme(line_colors = c("darkblue", "darkred")),
                plot_title = "US Monthly Tourist Arrivals & Spending, seasonally adjusted (SAAR) (2000-2019)")

2. Time Series Analysis

a. Assumption and Validation

To conduct a rigorous statistical analysis of U.S. monthly tourist arrivals and spending, we assume that each account in our data is independent of one another and meets all criteria described by the International Trade Administration for their corresponding reporting data. The following paragraphs outline our general approach to finding the best models for our data.

• Before fitting any model, we will visualize the data with a time series plot and compute the Autocorrelation Function (ACF) and Partial Autocorrelation Function (PACF) plots to understand the underlying patterns if appropriate. We will then perform several stationarity tests, including the Augmented Dickey-Fuller (ADF) test, Phillips-Perron (PP) test, Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test, and the Elliott-Rothenberg-Stock (ERS) test. More information on these tests and beyond will be available in the Appendix.

• After fitting the models, we will validate their effectiveness by plotting the fitted values with the original data and then analyzing the residuals. We will check for independence, normality, and stationarity of the residuals, employing tests such as the Shapiro-Wilk, Anderson-Darling, Jarque-Bera, and Lilliefors tests to ensure they follow a normal distribution. Additionally, ACF and density plots, Box-Ljung test, and other stationarity tests will inform us about potential issues with correlation and stationarity. In general, if the residuals are stationary, independent, and normally distributed, the model’s predictions can be deemed reliable.

• When comparing models, we will utilize the Akaike Information Criterion (AIC), Bayesian Information Criterion (BIC), and R-squared values (for linear models). A lower AIC and BIC, along with a higher R-squared value, indicate a better-fitting model. If residuals are stationary, we may enhance the model by incorporating additional terms to improve its predictive capability.

• For validation, we will implement a cross-validation approach by excluding data from 2019, training our models with the available data up to 2018, and subsequently comparing the model’s predictions against actual arrivals and spending for that year. Finally, we will forecast tourist arrivals and spending for two years into the future, assuming the absence of COVID-19-related disruptions.

b. Function of Time Models

In the event of identifying non-stationarity due to trends within the time series, we will apply function of time models to detrend the data. Specifically, we will utilize regression techniques to fit linear, quadratic, sinusoidal, and periodic trends with seasonal means. By identifying the appropriate model that best captures these trends, we will aim to derive residuals that should exhibit stationary behavior—constant variance and no serial correlation. Within this regard, we will also utilize decomposition techniques and plot accordingly to compare with our function of time models.

c. SARIMA Models

Seasonal Autoregressive Integrated Moving Average (SARIMA) models will be employed to account for seasonal behavior in the data. SARIMA models are denoted as SARIMA(p, d, q) × (P, D, Q) [s], where:

• p: Order of the autoregressive component

• d: Degree of differencing

• q: Order of the moving average component

• P: Seasonal autoregressive order

• D: Seasonal differencing order

• Q: Seasonal moving average order

• s: Length of the seasonal cycle

Prior to fitting a SARIMA model, we will first difference the data (normally and/or seasonally) as necessary to ensure stationarity. Different models will be fitted (based on ACF and PACF plots), validated by series of statistical tests, and analyzed in terms of residuals to confirm that they are stationary and uncorrelated. If the residuals follow a normal distribution, we can confidently utilize the model for forecasting future values.

d. Cross Correlation Analysis

To investigate the relationship between exports and arrivals, we use cross-correlation (CCF). We will first attempt to find the best ARIMA fit for the exports time series, then filter the seasonally adjusted time series with coefficients found in the model. Lastly, we will carry out the ccf() function, visualizing potential lags in the relationship.

e. Special Functions

To facilitate this analysis, I have self-defined several custom functions in R to streamline various processes, including:

• check_stationarity(): This function performs stationarity tests (ADF, PP, KPSS, ERS) and summarizes the results in a user-friendly table.

  # Stationarity Tests
  check_stationarity <- function(data) {
    adf_result <- adf.test(data)
    pp_result <- pp.test(data)
    kpss_result <- kpss.test(data)
    urers_result <- ur.ers(data, type = "DF-GLS", model = "constant")
    results <- data.frame(
      Test = c("ADF Test", "Phillips-Perron Test", "KPSS Test", "ERS Test"),
      `Test Statistic` = c(round(adf_result[["statistic"]][["Dickey-Fuller"]], 3),
                           round(pp_result[["statistic"]][["Dickey-Fuller Z(alpha)"]], 3),
                           round(kpss_result[["statistic"]][["KPSS Level"]], 3),
                           round(urers_result@teststat, 3)),
      `P-Value` = c(round(adf_result[["p.value"]],2),
                    round(pp_result[["p.value"]],2),
                    round(kpss_result[["p.value"]],2),
                    "N/A"),
      Interpretation = c(
        ifelse(adf_result[["p.value"]] < 0.05, "Stationary", "Non-Stationary"),
        ifelse(pp_result[["p.value"]] < 0.05, "Stationary", "Non-Stationary"),
        ifelse(kpss_result[["p.value"]] < 0.05, "Non-Stationary", "Stationary"),
        ifelse(urers_result@teststat < urers_result@cval[2], "Stationary", "Non-Stationary")))
    kable(results, align = 'c') %>%
      kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed")) %>%
      add_header_above(c("Stationarity Test Summary" = 4))
  }

• four_plots(), two_plots(), and diff_plots(): These functions generate ACF and PACF plots for different segments of the data to visually assess autocorrelation.

  # ACF/PACF/First-Last Third Plotting
  four_plots <- function(data) {
    n <- length(data)
    first_third <- data[1:floor(n / 3)]
    last_third <- data[floor(2 * n / 3):n]
    par(mfrow = c(2, 2))
    acf(data, main = "ACF of Whole Data", lag.max =30)
    pacf(data, main = "PACF of Whole Data", lag.max =30)
    acf(first_third, main = "ACF of First Third")
    acf(last_third, main = "ACF of Last Third")
    par(mfrow = c(1, 1))
  }
  
  # ACF/PACF Plotting
  two_plots <- function(data) {
    par(mfrow = c(1, 2))
    acf(data, main = "ACF of Whole Data", lag.max =50)
    pacf(data, main = "PACF of Whole Data", lag.max =50)
    par(mfrow = c(1, 1))
  }
  
  # Differencing Plotting
  diff_plots <- function(data, name) {
    layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
    astsa::tsplot(data, main= name)
    abline(h = mean(data, na.rm = TRUE), col = "red", lty = 2)
    acf(data, lag.max =30)
    pacf(data, lag.max =30)
    par(mfrow = c(1,1))
  }

• model_summary_table() and sarima_summary_table(): These functions summarize the key metrics of fitted models, including R^2 , AIC and BIC values, to facilitate comparison.

  # Model Summary
  model_summary_table <- function(models_list) {
    model_names <- names(models_list)
    if (is.null(model_names)) {
      model_names <- paste("Model", seq_along(models_list)) }
    summary_df <- data.frame(Model = character(),
                             `Adjusted  $R^2$ ` = numeric(),
                             AIC = numeric(),
                             BIC = numeric(),
                             stringsAsFactors = FALSE)
    for (i in seq_along(models_list)) {
      model <- models_list[[i]]
      model_summary <- summary(model)
      summary_df[i, ] <- c(
        model_names[i],
        round(model_summary[["adj.r.squared"]], 2),
        round(stats::AIC(model), 2),
        round(stats::BIC(model), 2)) }
    kable(summary_df, align = 'c', col.names = c("Model", "Adjusted R²", "AIC", "BIC")) %>%
      kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed")) %>%
      add_header_above(c("Model Summary" = 4))
  }
  
  
  # SARIMA Model Summary
  sarima_summary_table <- function(models_list) {
    model_names <- names(models_list)
    if (is.null(model_names)) {
      model_names <- paste("Model", seq_along(models_list)) }
    summary_df <- data.frame(Model = character(),
                             AIC = numeric(),
                             BIC = numeric(),
                             stringsAsFactors = FALSE)
    for (i in seq_along(models_list)) {
      model <- models_list[[i]]
      model_summary <- summary(model)
      summary_df[i, ] <- c(
        model_names[i],
        round(stats::AIC(model), 2),
        round(stats::BIC(model), 2)) }
    kable(summary_df, align = 'c', col.names = c("Model",  "AIC", "BIC")) %>%
      kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed")) %>%
      add_header_above(c("Model Summary" = 3))
  }

• check_normality_residuals(): This function evaluates the residuals from fitted models through a series of different tests to ensure they conform to normality assumptions.

  # Normality Tests Function
  check_normality_residuals <- function(model) {
    residuals <- residuals(model)
    ad_result <- ad.test(residuals)
    jb_result <- jarque.bera.test(residuals)
    lillie_result <- lillie.test(residuals)
    skew <- skewness(residuals)
    kurt <- kurtosis(residuals)
    # Summary table
    results <- data.frame(
      Test = c("Anderson-Darling Test", "Jarque-Bera Test", "Lilliefors Test", "Skewness", "Kurtosis"),
      `Test Statistic` = c(round(ad_result$statistic, 3),
                           round(jb_result$statistic, 3),
                           round(lillie_result$statistic, 3),
                           round(skew, 2),
                           round(kurt, 2)),
      `P-Value` = c(round(ad_result$p.value, 3),
                    round(jb_result$p.value, 3),
                    round(lillie_result$p.value, 3),
                    "N/A", "N/A"),
      Interpretation = c(
        ifelse(ad_result$p.value > 0.05, "Normal", "Non-Normal"),
        ifelse(jb_result$p.value > 0.05, "Normal", "Non-Normal"),
        ifelse(lillie_result$p.value > 0.05, "Normal", "Non-Normal"),
        ifelse(abs(skew) < 0.5, "Approximately Symmetric", ifelse(skew > 0, "Right-Skewed", "Left-Skewed")),
        ifelse(abs(kurt - 3) < 0.5, "Mesokurtic", ifelse(kurt > 3, "Leptokurtic (Heavy Tails)", "Platykurtic (Light Tails)"))
      )
    )
    kable(results, align = 'c') %>%
      kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed")) %>%
      add_header_above(c("Normality Test Summary" = 4))
  }

• ts_summary_table(): This function summarizes key characteristics of time series data.

  # Time Series Summary Table
  ts_summary_table <- function(ts_list) {
  
    ts_names <- names(ts_list)
    if (is.null(ts_names)) {
      ts_names <- paste("Series", seq_along(ts_list))  }
  
    summary_df <- data.frame(Series = character(),
                             Start = character(),
                             End = character(),
                             Min = numeric(),
                             Max = numeric(),
                             Mean = numeric(),
                             Median = numeric(),
                             stringsAsFactors = FALSE)
  
    for (i in seq_along(ts_list)) {
      ts_data <- ts_list[[i]]
      summary_df[i, ] <- c(
        ts_names[i],
        start(ts_data) %>% paste(collapse = "-"),
        end(ts_data) %>% paste(collapse = "-"),
        round(min(ts_data, na.rm = TRUE), 2),
        round(max(ts_data, na.rm = TRUE), 2),
        round(mean(ts_data, na.rm = TRUE), 2),
        round(median(ts_data, na.rm = TRUE), 2)
      )
    }
  
    kable(summary_df, align = 'c', 
          col.names = c("Series", "Start", "End", "Min", "Max", "Mean", "Median")) %>%
      kable_styling(full_width = F, bootstrap_options = c("striped", "hover", "condensed")) %>%
      add_header_above(c("Time Series Summary" = 7))
  }

• calculate_r2(): This function calculates the R^2 value for model predictions against actual values.

  calculate_r2 <- function(predicted_values, actual_values) {
    ss_residual <- sum((actual_values - predicted_values)^2)
    ss_total <- sum((actual_values - mean(actual_values))^2)
    r_squared <- 1 - (ss_residual / ss_total)
  
    cat("R-squared of the model is:", round(r_squared, 4), "\n")}

These functions enhance the analysis process by providing clarity and efficiency. Each function’s details and the tests used are elaborated upon in the appendix.

III. Results & Conclusions

1. Descriptive Analysis

a. Tourist Arrivals

U.S. monthly tourist arrivals showed a notable increase over the observed period, beginning with 2,866,229 visitors in January 2000 and reaching 6,899,661 by July 2024, indicating an overall upward trend. The plot below illustrates the monthly arrival data from 2000 to July 2024, segmented by different world regions as well as the total aggregate. The mean monthly arrivals were approximately 4,709,139, with a median of 4,817,910. Notably, Western Europe accounted for an average of 902,992 monthly arrivals, Asia for 618,172, and South America for 301,831.

The lowest recorded arrivals occurred in April 2020, when the number of visitors dropped sharply to 248,486, corresponding to the onset of U.S. border closures and social distancing measures due to the COVID-19 pandemic (Gonzalez-Barrera, 2020). Conversely, the highest number of monthly arrivals was recorded in August 2014, with 8,418,370 visitors.

astsa::tsplot(x= Masterdata0$Time, y= Masterdata0$Total/1000000,
       col = "darkblue",
       lwd = 2,
       main = "US Monthly Tourist Arrivals (2000-2024)",
       xlab = "Year",
       ylab = "Tourist Arrivals (millions)")
  
  lines(x= Masterdata0$Time, y= Masterdata0$Africa/1000000, col = 1, lwd = 1.5)
  lines(x= Masterdata0$Time, y= Masterdata0$Asia/1000000, col = 2, lwd = 1.5)
  lines(x= Masterdata0$Time, y= Masterdata0$W.Europe/1000000, col = 3, lwd = 1.5)
  lines(x= Masterdata0$Time, y= Masterdata0$E.Europe/1000000, col = 4, lwd = 1.5)
  lines(x= Masterdata0$Time, y= Masterdata0$C.America/1000000, col = 5, lwd = 1.5)
  lines(x= Masterdata0$Time, y= Masterdata0$Caribbean/1000000, col = 6, lwd = 1.5)
  lines(x= Masterdata0$Time, y= Masterdata0$S.America/1000000, col = 7, lwd = 1.5)
  lines(x= Masterdata0$Time, y= Masterdata0$Oceania/1000000, col = 8, lwd = 1.5)
  lines(x= Masterdata0$Time, y= Masterdata0$M.East/1000000, col = 9, lwd = 1.5)
  
  legend("topleft", col=c("darkblue",1:9), lty=1, lwd=2, pch=20, 
         legend=c("Total", "Africa", "Asia", "W.Europe", "E.Europe", "C.America", "Caribbean", "S.America", "Oceania","M.East"), bg="white",ncol=2)

Seasonal patterns are evident in the plot above, with clear peaks and troughs repeating annually. The plot below signifies that the summer months of July and August consistently show the highest arrivals, while January and February experience the lowest. This trend aligns with travel behavior, as summer weather is more favorable for tourism, and winter months see a decrease in travel as people return to work or school after the holidays. Additionally, colder weather and snow in northern regions may deter travel during these months.

p <- ggplot(Masterdata, aes(x=Month, y=Total/1000000, color = Year)) + 
  geom_line(alpha = 0.7) +
  labs(x = "Month",
       y = "Total Arrivals (millions)",
       title = "Total Arrivals by Month") +
  theme_minimal() + 
  scale_x_continuous(breaks=c(2,4, 6, 8,10,12)) +
  transition_states(Year) 

animate(p, renderer = gifski_renderer())

b. Tourist Spending

Monthly tourist spending, or “Tourist Exports,” also shows a marked increase, starting at 9.2 billion dollars in January 2000 and reaching 21.1 billion dollars in July 2024. The mean monthly spending over this period was 13.7 billion dollars, with a median of 13.1 billion dollars. The lowest recorded spending occurred in September 2020, at 3.7 billion dollars, while the highest was in May 2024, with 21.66 billion dollars. However, it is noteworthy that this dataset is seasonally adjusted, which means that nominal spending figures may vary; for instance, the lowest nominal spending month may have occurred earlier in 2020 when the tourism sector was effectively shut down.

ggplot(Masterdata0, aes(x = Time, y = Exports)) +
  geom_line(color = "darkred", size = 1) +
  ggtitle("US Monthly Tourist Exports, seasonally adjusted (2000-2024)") +
  xlab("Month") +
  ylab("Export Values (million $)") +
  scale_x_date(date_labels = "%b %Y", date_breaks = "36 months")+
  theme_light()+
  scale_y_continuous(labels = scales::label_number(scale_cut = cut_short_scale()))

c. Time Series summary

The table below compares statistical summaries of monthly tourist arrivals and exports for the entire period (2000-2024) versus the period excluding data from January 2020 onward.

  ar_full <- ts(Masterdata0$Total,start= c(2000,1), frequency = 12)
  ex_full <- ts(Masterdata0$Exports,start= c(2000,1), frequency = 12)
  ts_series <- list("Tourist Arrivals (2000-2024) (millions)" = ar_full/1000000,
                    "Tourist Arrivals (2000-2019) (millions)" = Arrivals_ts/1000000,
                    "Tourist Spending (2000-2024) (billions USD)" = ex_full/1000,
                    "Tourist Spending (2000-2019) (billions USD)" = Exports_ts/1000)

ts_summary_table(ts_series)

Time Series Summary
Series	Start	End	Min	Max	Mean	Median
Tourist Arrivals (2000-2024) (millions)	2000-1	2024-7	0.25	8.42	4.71	4.82
Tourist Arrivals (2000-2019) (millions)	2000-1	2019-12	2.09	8.42	4.96	4.91
Tourist Spending (2000-2024) (billions USD)	2000-1	2024-7	3.72	21.66	13.7	13.11
Tourist Spending (2000-2019) (billions USD)	2000-1	2019-12	6.7	20.82	13.86	12.89

2. Function of Time Models

a. Tourist Arrivals

In total, we evaluated six different models (four of which will be shown) to capture trends in the total arrivals time series, incorporating linear, quadratic, seasonal, sinusoidal components, and lagged values. In order to mitigate the multicollinearity problems associated with VIF, we intentionally tried to use the centered time object in all of our models, if appropriate.

Linear

This model shows a clear upward trend in the observed data. However, in terms of fitting, with dependent residuals, we cannot trust any metrics, like R^2, AIC, or BIC, nor the forecasts. Non-normality also affects the reliability of p-values and confidence intervals. Not a good fit, overall.

After trying this model, we attempted to add a quadratic term into the model, which yields not statistically significant coefficient. Hence, all subsequent models will not include any quadratic term.

t <- time(Arrivals_ts)
t_centered <- t - mean(t)

# Linear
ar_lm <- lm(Arrivals_ts ~ t_centered)
summary(ar_lm)

Call:
lm(formula = Arrivals_ts ~ t_centered)

Residuals:
     Min       1Q   Median       3Q      Max 
-1828757  -525528   -55961   354770  2465742 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4958137      50483   98.22   <2e-16 ***
t_centered    215025       8744   24.59   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 782100 on 238 degrees of freedom
Multiple R-squared:  0.7176,    Adjusted R-squared:  0.7164 
F-statistic: 604.7 on 1 and 238 DF,  p-value: < 2.2e-16

ts.plot(cbind(Arrivals_ts/1000000,ar_lm$fitted.values/1000000), col=c("darkblue", "darkred"),
       ylab = "Arrivals (millions)", main = "Original Data with Fitted Values - Linear", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"), col=c("darkblue", "darkred"), lty=1:1,lwd=2)

checkresiduals(ar_lm$residuals/1000000)

    Ljung-Box test

data:  Residuals
Q* = 176.1, df = 10, p-value < 2.2e-16

Model df: 0.   Total lags used: 10

check_stationarity(ar_lm$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-5.456	0.01	Stationary
Phillips-Perron Test	-113.997	0.01	Stationary
KPSS Test	0.109	0.1	Stationary
ERS Test	-7.740	N/A	Stationary

check_normality_residuals(ar_lm)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	2.066	0	Non-Normal
Jarque-Bera Test	13.541	0.001	Non-Normal
Lilliefors Test	0.084	0	Non-Normal
Skewness	0.570	N/A	Right-Skewed
Kurtosis	3.250	N/A	Mesokurtic

Linear + Seasonal

This model has fitted values matching observed data quite well. However, as the residuals lack independence, stationarity, randomness, and normality, we cannot trust R^2, AIC, or BIC values, nor p-values, confidence intervals, and the forecasts. A better, but not a good fit, yet.

# Linear + Seasonal Trend 
t_month <- factor(month(Masterdata$Time))
ar_month <- lm(Arrivals_ts ~ t_centered + t_month)
summary(ar_month)

Call:
lm(formula = Arrivals_ts ~ t_centered + t_month)

Residuals:
     Min       1Q   Median       3Q      Max 
-1123838  -251167    14405   276126  1177669 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4279501      96805  44.207  < 2e-16 ***
t_centered    213315       4845  44.027  < 2e-16 ***
t_month2     -402561     136868  -2.941  0.00361 ** 
t_month3      594818     136870   4.346 2.09e-05 ***
t_month4      845284     136873   6.176 3.01e-09 ***
t_month5      710543     136877   5.191 4.64e-07 ***
t_month6      451356     136882   3.297  0.00113 ** 
t_month7     1854503     136889  13.547  < 2e-16 ***
t_month8     2018620     136897  14.746  < 2e-16 ***
t_month9      678724     136906   4.958 1.40e-06 ***
t_month10     622646     136916   4.548 8.83e-06 ***
t_month11     181691     136927   1.327  0.18587    
t_month12     587998     136940   4.294 2.60e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 432800 on 227 degrees of freedom
Multiple R-squared:  0.9175,    Adjusted R-squared:  0.9131 
F-statistic: 210.4 on 12 and 227 DF,  p-value: < 2.2e-16

ts.plot(cbind(Arrivals_ts/1000000,ar_month$fitted.values/1000000 ), col=c("darkblue", "darkred"),
       ylab = "Arrivals (millions)", main = "Original Data with Fitted Values - Linear + Seasonal", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"), col=c("darkblue", "darkred"), lty=1:1,lwd=2)

checkresiduals(ar_month$residuals/1000000)

    Ljung-Box test

data:  Residuals
Q* = 713.33, df = 10, p-value < 2.2e-16

Model df: 0.   Total lags used: 10

check_stationarity(ar_month$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-3.270	0.08	Non-Stationary
Phillips-Perron Test	-44.256	0.01	Stationary
KPSS Test	0.233	0.1	Stationary
ERS Test	-1.551	N/A	Non-Stationary

check_normality_residuals(ar_month)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	0.757	0.048	Non-Normal
Jarque-Bera Test	0.225	0.894	Normal
Lilliefors Test	0.063	0.022	Non-Normal
Skewness	0.010	N/A	Approximately Symmetric
Kurtosis	3.150	N/A	Mesokurtic

Linear + Seasonal + Cosine

This model alsoo has fitted values matching observed data quite well. However, as the residuals lack independence and stationarity, we cannot trust R^2, AIC, or BIC values, nor p-values, confidence intervals, and the forecasts. A slightly better, but not a great fit, yet.

# Linear + Seasonal + Cosine

cosine_t <- cosine_t <- cos(2 * pi * t / 12)

ar_both <- lm(Arrivals_ts ~ t_centered  + cosine_t + t_month)
summary(ar_both)

Call:
lm(formula = Arrivals_ts ~ t_centered + cosine_t + t_month)

Residuals:
     Min       1Q   Median       3Q      Max 
-1035888  -297866     9284   274209  1196632 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4297509      95791  44.864  < 2e-16 ***
t_centered    213407       4782  44.625  < 2e-16 ***
cosine_t     -111170      41988  -2.648 0.008676 ** 
t_month2     -402376     135092  -2.979 0.003212 ** 
t_month3      595153     135093   4.405 1.63e-05 ***
t_month4      845736     135096   6.260 1.91e-09 ***
t_month5      711075     135100   5.263 3.28e-07 ***
t_month6      451933     135106   3.345 0.000963 ***
t_month7     1855090     135112  13.730  < 2e-16 ***
t_month8     2019183     135120  14.944  < 2e-16 ***
t_month9      679226     135129   5.027 1.02e-06 ***
t_month10     623052     135139   4.610 6.72e-06 ***
t_month11     181965     135150   1.346 0.179523    
t_month12     588107     135162   4.351 2.05e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 427200 on 226 degrees of freedom
Multiple R-squared:   0.92, Adjusted R-squared:  0.9154 
F-statistic: 199.9 on 13 and 226 DF,  p-value: < 2.2e-16

ts.plot(cbind(Arrivals_ts/1000000,ar_both$fitted.values/1000000 ), col=c("darkblue", "darkred"),
       ylab = "Arrivals (millions)", main = "Original Data with Fitted Values - Linear + Seasonal + Cosine", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"), col=c("darkblue", "darkred"), lty=1:1,lwd=2)

checkresiduals(ar_both$residuals)

    Ljung-Box test

data:  Residuals
Q* = 694.36, df = 10, p-value < 2.2e-16

Model df: 0.   Total lags used: 10

check_stationarity(ar_both$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-3.193	0.09	Non-Stationary
Phillips-Perron Test	-45.899	0.01	Stationary
KPSS Test	0.226	0.1	Stationary
ERS Test	-1.676	N/A	Non-Stationary

check_normality_residuals(ar_both)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	0.229	0.809	Normal
Jarque-Bera Test	0.765	0.682	Normal
Lilliefors Test	0.035	0.665	Normal
Skewness	0.130	N/A	Approximately Symmetric
Kurtosis	2.890	N/A	Mesokurtic

Linear + Seasonal + Lag 1

This model fits the data the best, even too good to feel true, as R^2, both normal and adjusted, are 1, which raises concerns about over-fitting. In terms of residuals, though not-normal (highly due to the presence of outliers), they are independent (passibng the L-jung Box Test, though might still have significant lags of 11 and 12) and stationary. As the model might be over-fitting, we will not trust this model much.

# Linear + Seasonal + Lag 1

ar_1 <- stats::lag(Arrivals_ts, 1)
ar_final <-lm(Arrivals_ts ~ t_centered + t_month + ar_1)
summary(ar_final)

Call:
lm(formula = Arrivals_ts ~ t_centered + t_month + ar_1)

Residuals:
       Min         1Q     Median         3Q        Max 
-5.116e-09 -3.440e-11  3.100e-12  4.740e-11  1.037e-09 

Coefficients:
              Estimate Std. Error    t value Pr(>|t|)    
(Intercept)  3.801e-09  2.527e-10  1.504e+01  < 2e-16 ***
t_centered   1.895e-10  1.260e-11  1.504e+01  < 2e-16 ***
t_month2    -3.575e-10  1.174e-10 -3.045e+00 0.002605 ** 
t_month3     5.283e-10  1.200e-10  4.404e+00 1.64e-05 ***
t_month4     7.508e-10  1.246e-10  6.027e+00 6.73e-09 ***
t_month5     6.311e-10  1.219e-10  5.176e+00 4.99e-07 ***
t_month6     4.009e-10  1.180e-10  3.397e+00 0.000804 ***
t_month7     1.647e-09  1.550e-10  1.062e+01  < 2e-16 ***
t_month8     1.793e-09  1.613e-10  1.112e+01  < 2e-16 ***
t_month9     6.028e-10  1.214e-10  4.967e+00 1.34e-06 ***
t_month10    5.530e-10  1.204e-10  4.592e+00 7.29e-06 ***
t_month11    1.614e-10  1.158e-10  1.394e+00 0.164650    
t_month12    5.222e-10  1.199e-10  4.355e+00 2.02e-05 ***
ar_1         1.000e+00  5.589e-17  1.789e+16  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.645e-10 on 226 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 2.985e+32 on 13 and 226 DF,  p-value: < 2.2e-16

ts.plot(cbind(Arrivals_ts/1000000,ar_final$fitted.values/1000000 ), col=c("darkblue", "darkred"), 
       ylab = "Arrivals (millions)", main = "Original Data with Fitted Values - Linear + Seasonal + Lag1", lwd =c(2.5,1))
legend("topleft", legend=c("Original Data", "Fitted Values"), col=c("darkblue", "darkred"), lty=1:1,lwd=2)

checkresiduals(ar_final$residuals)

    Ljung-Box test

data:  Residuals
Q* = 1.2514, df = 10, p-value = 0.9995

Model df: 0.   Total lags used: 10

check_stationarity(ar_final$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-6.353	0.01	Stationary
Phillips-Perron Test	-234.429	0.01	Stationary
KPSS Test	0.107	0.1	Stationary
ERS Test	-2.244	N/A	Stationary

check_normality_residuals(ar_final)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	49.482	0	Non-Normal
Jarque-Bera Test	329726.796	0	Non-Normal
Lilliefors Test	0.358	0	Non-Normal
Skewness	-12.520	N/A	Left-Skewed
Kurtosis	182.850	N/A	Leptokurtic (Heavy Tails)

Model Comparison

Of all of the models we fit, based on stationarity and independece, we can only trust the R^2, AIC, and BIC values of the last model with the lagged term. Yet, we need to proceed with caution as over-fitting can be a valid concern.

# Summary
ar_models <- list("Linear Model" = ar_lm,
                  "Linear + Seasonal" = ar_month,
                  "Linear + Seasonal + Cosine" = ar_both,
                  "Linear + Seasonal + Lag 1" = ar_final)
model_summary_table(ar_models)

Model Summary
Model	Adjusted R²	AIC	BIC
Linear Model	0.72	7198.54	7208.98
Linear + Seasonal	0.91	6925.19	6973.92
Linear + Seasonal + Cosine	0.92	6919.86	6972.07
Linear + Seasonal + Lag 1	1	-9734.97	-9682.76

Decomposition

Here, we try to decompose the time series. As shown in the plots, there are significant seasonality of frequency 12 months.

# Classical Seasonal Decomposition by Moving Averages
plot(decompose(Arrivals_ts/1000000))

# STL Decomposition
plot(stl(Arrivals_ts/1000000, s.window = 12))

Cross-Validation

Here, we trained the model with data up to 2018, and then tested against actual data in 2019. This technique is called Leave One Out Cross Validation (LOOCV) (Olivo, 2024). The predicted values did not capture much variablity in 2019 and the \(R^2\) is negative, suggesting a poor fit.

# Train on Data Up to 2018 and Predict for 2019 (with 95% CI)
train_data <- window(Arrivals_ts, end = c(2018, 12))
test_data <- window(Arrivals_ts, start = c(2019, 1), end = c(2019, 12))
t_train <- time(train_data)
t_month_train <- t_month[1:228]
cosine_t_train <- cosine_t[1:228]
model_train <- lm(train_data ~ t_train + cosine_t_train + t_month_train)
t_test <- time(test_data)
t_month_test <- t_month[229:240]
cosine_t_test <- cosine_t[229:240]
predictions_2019 <- predict(model_train,
                            newdata = data.frame(t_train = t_test,
                                                 cosine_t_train = cosine_t_test,
                                                 t_month_train = t_month_test),
                            interval = "confidence", level = 0.95)

# Combine actual and predicted data
time_full_2019 <- data.frame(
  Date = c(time(train_data), time(test_data)),
  Actual = c(as.numeric(train_data), as.numeric(test_data)),
  Predicted = c(rep(NA, length(train_data)), predictions_2019[,"fit"]),
  Lower_CI = c(rep(NA, length(train_data)), predictions_2019[,"lwr"]),
  Upper_CI = c(rep(NA, length(train_data)), predictions_2019[,"upr"]))

time_only_2019 <- data.frame(
  Date = as.Date(time(test_data)),
  Actual = test_data,
  Predicted = predictions_2019[,"fit"],
  Lower_CI = predictions_2019[,"lwr"],
  Upper_CI = predictions_2019[,"upr"])

ar_pred_2019 <- ggplot(time_only_2019, aes(x = Date)) +
  geom_line(aes(y = Actual/1000000, color = "Actual")) +
  geom_line(aes(y = Predicted/1000000, color = "Predicted")) +
  geom_ribbon(aes(ymin = Lower_CI/1000000, ymax = Upper_CI/1000000), fill = "grey30", alpha = 0.4) +
  labs(title = "Actual and Predicted Arrivals (2019) with 95% CI",
       y = "Tourist Arrivals (millions)")+
  scale_x_date(date_labels = "%b %Y", date_breaks = "2 months") +
  theme_minimal() +
  scale_color_manual(values = c("Actual" = "darkblue", "Predicted" = "red"))

ar_pred_full <- ggplot(time_full_2019, aes(x = Date)) +
  geom_line(aes(y = Actual/1000000, color = "Actual")) +
  geom_line(aes(y = Predicted/1000000, color = "Predicted")) +
  geom_ribbon(aes(ymin = Lower_CI/1000000, ymax = Upper_CI/1000000), fill = "grey30", alpha = 0.4) +
  labs(title = "Actual and Predicted Arrivals (2019) with 95% CI",
       y = "Tourist Arrivals (millions)") +
  theme_minimal() +
  scale_color_manual(values = c("Actual" = "darkblue", "Predicted" = "red"))
grid.arrange(ar_pred_full, ar_pred_2019, nrow=2)

# R^2
calculate_r2(as.numeric(predictions_2019),as.numeric(test_data))

R-squared of the model is: -0.675 

Prediction

Based on the poor result in the cross-validation test, when we applied the model to forecast arrivals for hypothetical years 2020 and 2021 (assuming no impact from COVID-19), the predictions would likely not help much.

# Fit Model on Full Data and Forecast Two Years Ahead (with 95% CI)
ar_pred <- lm(Arrivals_ts ~ t  + cosine_t + t_month)
future_dates <- ts(rep(NA, 24),start= c(2020,1), frequency = 12)
t_future <- time(future_dates)
cosine_t_future <- cos(2 * pi * time(future_dates) / 12)
t_month_future <- factor(month(seq(as.Date("2020-01-01"),
                                   as.Date("2021-12-01"), by = "month")))
predictions_future <- predict(ar_pred,
                              newdata = data.frame(t = t_future,
                                                   cosine_t = cosine_t_future,
                                                   t_month = t_month_future),
                              interval = "confidence", level = 0.95)

# Combine actual data with forecast
forecast_2020_2021 <- data.frame(
  Date = c(time(Arrivals_ts), time(future_dates)),
  Actual = c(as.numeric(Arrivals_ts), rep(NA, length(future_dates))),
  Forecast = c(rep(NA, length(Arrivals_ts)), predictions_future[,"fit"]),
  Lower_CI = c(rep(NA, length(Arrivals_ts)), predictions_future[,"lwr"]),
  Upper_CI = c(rep(NA, length(Arrivals_ts)), predictions_future[,"upr"]))

# Plot
ggplot(forecast_2020_2021, aes(x = Date)) +
  geom_line(aes(y = Actual/1000000, color = "Actual")) +
  geom_line(aes(y = Forecast/1000000, color = "Forecast")) +
  geom_ribbon(aes(ymin = Lower_CI/1000000, ymax = Upper_CI/1000000), fill = "grey30", alpha = 0.4) +
  labs(title = "Tourist Arrivals Forecast (2020-2021) with 95% CI",
       y = "Tourist Arrivals (millions)") +
  theme_minimal() +
  scale_color_manual(values = c("Actual" = "darkblue", "Forecast" = "red"))

b. Tourist Spending

Linear

This model exhibits a clear upward trend over time. However, in terms of fitting, with dependent, non-stationary and non-normal residuals, we cannot trust any metrics, like R^2, AIC, or BIC, nor the forecasts. Non-normality also affects the reliability of p-values and confidence intervals. Not a good fit, overall.

t <- time(Exports_ts)
t_centered <- t - mean(t)


# Linear
ex_lm <- lm(Exports_ts ~ t_centered)
summary(ex_lm)

Call:
lm(formula = Exports_ts ~ t_centered)

Residuals:
    Min      1Q  Median      3Q     Max 
-2551.6  -929.9  -169.9   701.6  3326.2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 13864.95      80.50  172.23   <2e-16 ***
t_centered    736.80      13.94   52.84   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1247 on 238 degrees of freedom
Multiple R-squared:  0.9215,    Adjusted R-squared:  0.9211 
F-statistic:  2792 on 1 and 238 DF,  p-value: < 2.2e-16

ts.plot(cbind(Exports_ts/1000,ex_lm$fitted.values/1000), 
        col=c("darkblue", "darkred"), ylab = "Exports (billions $)", 
        main = "Original Data with Fitted Values - Linear", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"), 
       col=c("darkblue", "darkred"), lty=1:1,lwd=2)

checkresiduals(ex_lm$residuals)

    Ljung-Box test

data:  Residuals
Q* = 1464.4, df = 10, p-value < 2.2e-16

Model df: 0.   Total lags used: 10

check_stationarity(ex_lm$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-3.005	0.15	Non-Stationary
Phillips-Perron Test	-8.956	0.61	Non-Stationary
KPSS Test	0.515	0.04	Non-Stationary
ERS Test	-0.496	N/A	Non-Stationary

check_normality_residuals(ex_lm)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	2.712	0	Non-Normal
Jarque-Bera Test	15.144	0.001	Non-Normal
Lilliefors Test	0.085	0	Non-Normal
Skewness	0.620	N/A	Right-Skewed
Kurtosis	3.010	N/A	Mesokurtic

Quadratic

This model appears to be slightly better, but all of the problems with residuals remain, suggesting we cannot trust much of the model’s predictive power.

# Quadratic
ex_quad <- lm(Exports_ts ~ t_centered + I(t^2) )
summary(ex_quad)

Call:
lm(formula = Exports_ts ~ t_centered + I(t^2))

Residuals:
    Min      1Q  Median      3Q     Max 
-2504.8  -837.3  -230.1   968.2  2436.5 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -5.899e+07  1.024e+07  -5.762 2.57e-08 ***
t_centered  -5.797e+04  1.019e+04  -5.691 3.71e-08 ***
I(t^2)       1.461e+01  2.534e+00   5.763 2.55e-08 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1170 on 237 degrees of freedom
Multiple R-squared:  0.9311,    Adjusted R-squared:  0.9305 
F-statistic:  1602 on 2 and 237 DF,  p-value: < 2.2e-16

ts.plot(cbind(Exports_ts/1000,ex_quad$fitted.values/1000), 
        col=c("darkblue", "darkred"), ylab = "Exports (billions $)", 
        main = "Original Data with Fitted Values - Quadratic", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"), 
       col=c("darkblue", "darkred"), lty=1:1,lwd=2)

checkresiduals(ex_quad$residuals)

    Ljung-Box test

data:  Residuals
Q* = 1407.8, df = 10, p-value < 2.2e-16

Model df: 0.   Total lags used: 10

check_stationarity(ex_quad$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-2.229	0.48	Non-Stationary
Phillips-Perron Test	-7.588	0.68	Non-Stationary
KPSS Test	0.421	0.07	Stationary
ERS Test	-0.611	N/A	Non-Stationary

check_normality_residuals(ex_quad)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	2.727	0	Non-Normal
Jarque-Bera Test	8.759	0.013	Non-Normal
Lilliefors Test	0.098	0	Non-Normal
Skewness	0.220	N/A	Approximately Symmetric
Kurtosis	2.180	N/A	Platykurtic (Light Tails)

Quadratic + Lag 1

This model does solve problems of dependence and non-staionarity in the residuals, but raises tremendous concerns of over-fitting. Hence, we do not wish to proceed with further analysis.

ex_1 <- stats::lag(Exports_ts, 1)
ex_final <-lm(Exports_ts ~ t_centered + I(t^2) + ex_1 )
summary(ex_final)

Call:
lm(formula = Exports_ts ~ t_centered + I(t^2) + ex_1)

Residuals:
       Min         1Q     Median         3Q        Max 
-2.100e-12 -8.140e-14 -9.000e-15  5.030e-14  1.477e-11 

Coefficients:
              Estimate Std. Error    t value Pr(>|t|)    
(Intercept)  5.239e-08  9.303e-09  5.632e+00 5.05e-08 ***
t_centered   5.149e-11  9.243e-12  5.571e+00 6.89e-08 ***
I(t^2)      -1.297e-14  2.303e-15 -5.633e+00 5.02e-08 ***
ex_1         1.000e+00  5.528e-17  1.809e+16  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 9.961e-13 on 236 degrees of freedom
Multiple R-squared:      1, Adjusted R-squared:      1 
F-statistic: 1.583e+33 on 3 and 236 DF,  p-value: < 2.2e-16

ts.plot(cbind(Exports_ts/1000,ex_final$fitted.values/1000), 
        col=c("darkblue", "darkred"), ylab = "Exports (billions $)", 
        main = "Original Data with Fitted Values - Quadratic + Lag1", lwd =c(2,1))
legend("topleft", legend=c("Original Data", "Fitted Values"), 
       col=c("darkblue", "darkred"), lty=1:1,lwd=2)

checkresiduals(ex_final$residuals)

    Ljung-Box test

data:  Residuals
Q* = 8.5456, df = 10, p-value = 0.5757

Model df: 0.   Total lags used: 10

check_stationarity(ex_final$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-3.516	0.04	Stationary
Phillips-Perron Test	-237.781	0.01	Stationary
KPSS Test	0.037	0.1	Stationary
ERS Test	-0.419	N/A	Non-Stationary

check_normality_residuals(ex_final)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	61.830	0	Non-Normal
Jarque-Bera Test	429668.502	0	Non-Normal
Lilliefors Test	0.402	0	Non-Normal
Skewness	13.840	N/A	Right-Skewed
Kurtosis	208.430	N/A	Leptokurtic (Heavy Tails)

Model Comparison

Here, a quadratic model emerged as the best fit, with an R-squared of 0.93, an AIC of 4077.33, and a BIC of 4091.25. However, this model’s residuals are not independent, suggesting poor reliability of these reported metrics.

# Summary
ex_models <- list("Linear Model" = ex_lm, 
                  "Quadratic Model" = ex_quad, 
                  "Quadratic + Lag1" = ex_final)
model_summary_table(ex_models)

Model Summary
Model	Adjusted R²	AIC	BIC
Linear Model	0.92	4106.8	4117.25
Quadratic Model	0.93	4077.33	4091.25
Quadratic + Lag1	1	-12577.73	-12560.32

Additionally, the quadratic model’s residuals, like those of the arrivals model, failed to meet stationarity, normality, and independence assumptions, limiting the reliability of its forecast intervals. It also possesses problem of collinearity.

check_model(ex_quad)

Decomposition

Here, we try to decompose the time series. As shown in the plots, while there seems to be signs of seasonality with frequency of 12 months, the magnitude of these seasonality signals are not so significant.

# Classical Seasonal Decomposition by Moving Averages
plot(decompose(Exports_ts/1000))

# STL Decomposition
plot(stl(Exports_ts/1000, s.window = 12))

Prediction

Based on all of the poor statistics and test results above with this quadratic model, we know for certain that this model’s predictive power should not be deemed reliable and all predictions have to be interpreted with caution. Hence, we skipped the cross-validation test and proceeded with the prediction.

# Fit Model on Full Data and Forecast Two Years Ahead (with 95% CI)
ex_pred <- lm(Arrivals_ts ~ t + I(t^2))


future_dates <- ts(rep(NA, 24), start = c(2020, 1), frequency = 12)
t_future <- time(future_dates)

predictions_future <- predict(ex_pred,
                              newdata = data.frame(t = t_future,
                                                   `I(t^2)` = t_future^2),
                              interval = "confidence", level = 0.95)


forecast_2020_2021 <- data.frame(
  Date = c(time(Arrivals_ts), time(future_dates)),
  Actual = c(as.numeric(Arrivals_ts), rep(NA, length(future_dates))),
  Forecast = c(rep(NA, length(Arrivals_ts)), predictions_future[,"fit"]),
  Lower_CI = c(rep(NA, length(Arrivals_ts)), predictions_future[,"lwr"]),
  Upper_CI = c(rep(NA, length(Arrivals_ts)), predictions_future[,"upr"]))


ggplot(forecast_2020_2021, aes(x = Date)) +
  geom_line(aes(y = Actual / 1000000, color = "Actual")) +
  geom_line(aes(y = Forecast / 1000000, color = "Forecast")) +
  geom_ribbon(aes(ymin = Lower_CI / 1000000, ymax = Upper_CI / 1000000), fill = "grey30", alpha = 0.4) +
  labs(title = "Tourist Arrivals Forecast (2020-2021) with 95% CI",
       y = "Tourist Arrivals (millions)") +
  theme_minimal() +
  scale_color_manual(values = c("Actual" = "darkblue", "Forecast" = "red"))

# Define Lagged Series
ex_1 <- stats::lag(Exports_ts, 1)

# Fit Model on Full Data Including Quadratic Trend and Lagged Term
ex_final <- lm(Exports_ts ~ t + I(t^2) + ex_1)

# Forecast Two Years Ahead (2020-2021) with 95% CI
future_dates <- ts(rep(NA, 24), start = c(2020, 1), frequency = 12)
t_future <- time(future_dates)

# Get the last observed lagged value for prediction purposes
last_observed_value <- tail(Exports_ts, 1)

# Prepare future lagged values by starting with the last observed value
lagged_future <- c(last_observed_value, rep(NA, length(t_future) - 1))

# Iterate to populate future lagged values based on predictions
for (i in 2:length(lagged_future)) {
  lagged_future[i] <- predict(ex_final,
                              newdata = data.frame(
                                t = t_future[i-1],
                                `I(t^2)` = t_future[i-1]^2,
                                ex_1 = lagged_future[i-1]
                              ))
}

# Generate future predictions
predictions_future <- predict(ex_final,
                              newdata = data.frame(
                                t = t_future,
                                `I(t^2)` = t_future^2,
                                ex_1 = lagged_future),
                              interval = "confidence", level = 0.95)

# Combine actual data with forecast
forecast_2020_2021 <- data.frame(
  Date = c(time(Exports_ts), time(future_dates)),
  Actual = c(as.numeric(Exports_ts), rep(NA, length(future_dates))),
  Forecast = c(rep(NA, length(Exports_ts)), predictions_future[,"fit"]),
  Lower_CI = c(rep(NA, length(Exports_ts)), predictions_future[,"lwr"]),
  Upper_CI = c(rep(NA, length(Exports_ts)), predictions_future[,"upr"])
)

# Plotting the Forecast
ggplot(forecast_2020_2021, aes(x = Date)) +
  geom_line(aes(y = Actual / 1000000, color = "Actual")) +
  geom_line(aes(y = Forecast / 1000000, color = "Forecast")) +
  geom_ribbon(aes(ymin = Lower_CI / 1000000, ymax = Upper_CI / 1000000), fill = "grey30", alpha = 0.4) +
  labs(title = "Exports Forecast (2020-2021) with 95% CI",
       y = "Exports (millions)") +
  theme_minimal() +
  scale_color_manual(values = c("Actual" = "darkblue", "Forecast" = "red"))

3. SARIMA Models

a. Tourist Arrivals

i. Achieving Stationarity

Time Series

This time series is definitely not stationary, based on the time series plot, the tests, and the ACF/PACF plots.

diff_plots(Arrivals_ts/1000000, "Tourist Arrivals, millions (2000-2019)")

check_stationarity(Arrivals_ts)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-5.456	0.01	Stationary
Phillips-Perron Test	-113.997	0.01	Stationary
KPSS Test	4.286	0.01	Non-Stationary
ERS Test	-1.360	N/A	Non-Stationary

Auto ARIMA

Auto.arima() suggests an ARIMA(3,0,2)(2,1,0)[12] with drift. From the outset, this model seems to fit the model well, with stationary and independent residuals. Yet, the residuals are not normal, suggesting poor reliability of forecasts. This step is only for reference, and we will attempt to build our only models.

# Auto.arima() - ARIMA(3,0,2)(2,1,0)[12] with drift
ar_initial <- auto.arima(Arrivals_ts)
ar_initial

Series: Arrivals_ts 
ARIMA(3,0,2)(2,1,0)[12] with drift 

Coefficients:
         ar1      ar2     ar3     ma1     ma2     sar1     sar2      drift
      0.8242  -0.2811  0.3483  -0.317  0.3689  -0.5835  -0.3680  13193.359
s.e.  0.2582   0.2985  0.1605   0.254  0.2163   0.0655   0.0682   6498.454

sigma^2 = 5.956e+10:  log likelihood = -3151.28
AIC=6320.56   AICc=6321.38   BIC=6351.42

ts.plot(cbind(Arrivals_ts/1000000,ar_initial$fitted/1000000), 
        col=c("darkblue", "darkred"), ylab = "Arrivals (millions)", 
        main = "Original Data with Fitted Values - Auto ARIMA", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"), 
       col=c("darkblue", "darkred"), lty=1:1,lwd=2)

checkresiduals(ar_initial)

    Ljung-Box test

data:  Residuals from ARIMA(3,0,2)(2,1,0)[12] with drift
Q* = 18.224, df = 17, p-value = 0.3748

Model df: 7.   Total lags used: 24

two_plots(ar_initial$residuals)

check_stationarity(ar_initial$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-5.533	0.01	Stationary
Phillips-Perron Test	-237.742	0.01	Stationary
KPSS Test	0.135	0.1	Stationary
ERS Test	-7.112	N/A	Stationary

check_normality_residuals(ar_initial)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	2.200	0	Non-Normal
Jarque-Bera Test	51.371	0	Non-Normal
Lilliefors Test	0.074	0.003	Non-Normal
Skewness	-0.230	N/A	Approximately Symmetric
Kurtosis	5.220	N/A	Leptokurtic (Heavy Tails)

Differencing

Recall the ACF plot with a slow decay in both seasonal and normal lags, we will try both techniques here.

# Seasonal Differencing
sd_ar <- diff(Arrivals_ts,12)
diff_plots(sd_ar/1000000,"Seasonally Differenced Time Series")

check_stationarity(sd_ar)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-3.026	0.14	Non-Stationary
Phillips-Perron Test	-70.739	0.01	Stationary
KPSS Test	0.305	0.1	Stationary
ERS Test	-4.092	N/A	Stationary

By only seasonally differencing, the series does not pass all stationarity tests.

# Normal Differencing
diff_ar <- diff(Arrivals_ts)
diff_plots(diff_ar/1000000, "Normally Differenced Time Series")

check_stationarity(diff_ar)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-12.989	0.01	Stationary
Phillips-Perron Test	-221.516	0.01	Stationary
KPSS Test	0.009	0.1	Stationary
ERS Test	-5.811	N/A	Stationary

By only normally differencing, the series is already stationary, based on the “by-eye” test and other statistical tests. Looking at the seasonal lags 12 and 24 decaying slowly in the ACF and the signicant lag 12 in the PACF, if we start fitting with only normal differencing, we will start with SAR 1. Note that the mean of residuals is zero, so we don’t need to include a drift.

# Normal & Seasonal Differencing
ar_nsd <- diff(diff(Arrivals_ts,12))
diff_plots(ar_nsd/1000000, "Normally & Seasonally Differenced Time Series")

check_stationarity(ar_nsd)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-6.957	0.01	Stationary
Phillips-Perron Test	-302.729	0.01	Stationary
KPSS Test	0.018	0.1	Stationary
ERS Test	-2.911	N/A	Stationary

By differencing both ways, we achieved a stationary series, with smaller variance and less autocorrelations. Based on the ACF and PACF, a next move could be MA 1. Note that the mean of residuals is zero, so we don’t need to include a drift.

ii. Model Fitting

SARIMA Attempt 1 (d=1, D=0)

In this attempt, we will only focus on SARIMA models with \(d=1, D=0\).

ARIMA 010100

First, we fit an ARIMA(0,1,0)(1,0,0)[12]. Based on the ACF and PACF, the best next move would be to add an AR1, or \(p=1\).

# ARIMA 010100
ar_010100 <- arima(Arrivals_ts, order=c(0,1,0),
                   seasonal=list(order=c(1,0,0), period=12))
ar_010100

Call:
arima(x = Arrivals_ts, order = c(0, 1, 0), seasonal = list(order = c(1, 0, 0), 
    period = 12))

Coefficients:
        sar1
      0.9001
s.e.  0.0238

sigma^2 estimated as 1.04e+11:  log likelihood = -3380.49,  aic = 6764.97

checkresiduals(ar_010100)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,0)(1,0,0)[12]
Q* = 109.87, df = 23, p-value = 2.663e-13

Model df: 1.   Total lags used: 24

two_plots(ar_010100$residuals)

check_stationarity(ar_010100$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-7.658	0.01	Stationary
Phillips-Perron Test	-310.831	0.01	Stationary
KPSS Test	0.015	0.1	Stationary
ERS Test	-7.790	N/A	Stationary

check_normality_residuals(ar_010100)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	1.812	0	Non-Normal
Jarque-Bera Test	34.803	0	Non-Normal
Lilliefors Test	0.069	0.008	Non-Normal
Skewness	-0.100	N/A	Approximately Symmetric
Kurtosis	4.850	N/A	Leptokurtic (Heavy Tails)

ARIMA 011100

Here, we fit an ARIMA(0,1,1)(1,0,0)[12]. Based on the ACF and PACF, the best next move would be to add an SMA1, or \(Q=1\).

# ARIMA 011100
ar_011100 <- arima(Arrivals_ts, order=c(0,1,1), 
                   seasonal=list(order=c(1,0,0), period=12))
ar_011100

Call:
arima(x = Arrivals_ts, order = c(0, 1, 1), seasonal = list(order = c(1, 0, 0), 
    period = 12))

Coefficients:
          ma1    sar1
      -0.4773  0.9175
s.e.   0.0585  0.0211

sigma^2 estimated as 8.338e+10:  log likelihood = -3355.34,  aic = 6716.69

checkresiduals(ar_011100)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,1)(1,0,0)[12]
Q* = 54.975, df = 22, p-value = 0.0001198

Model df: 2.   Total lags used: 24

two_plots(ar_011100$residuals)

check_stationarity(ar_011100$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-6.288	0.01	Stationary
Phillips-Perron Test	-251.548	0.01	Stationary
KPSS Test	0.022	0.1	Stationary
ERS Test	-7.745	N/A	Stationary

check_normality_residuals(ar_011100)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	1.573	0	Non-Normal
Jarque-Bera Test	40.733	0	Non-Normal
Lilliefors Test	0.067	0.011	Non-Normal
Skewness	0.040	N/A	Approximately Symmetric
Kurtosis	5.020	N/A	Leptokurtic (Heavy Tails)

ARIMA 011101

Here, we fit an ARIMA(0,1,1)(1,0,1)[12]. Based on the ACF, PACF, and other tests, we have finally found a model with independent and stationary residuals. All coefficients in the model are significant with small standard errors. However, the residuals are not normal; hence, we cannot trust the forecast intervals.

# ARIMA 011101
ar_011101 <- arima(Arrivals_ts, order=c(0,1,1), 
                   seasonal=list(order=c(1,0,1), period=12), method="ML")
ar_011101

Call:
arima(x = Arrivals_ts, order = c(0, 1, 1), seasonal = list(order = c(1, 0, 1), 
    period = 12), method = "ML")

Coefficients:
          ma1    sar1     sma1
      -0.4569  0.9929  -0.6560
s.e.   0.0605  0.0036   0.0523

sigma^2 estimated as 5.924e+10:  log likelihood = -3319.68,  aic = 6647.36

checkresiduals(ar_011101)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,1)(1,0,1)[12]
Q* = 21.741, df = 21, p-value = 0.4146

Model df: 3.   Total lags used: 24

two_plots(ar_011101$residuals)

check_stationarity(ar_011101$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-6.653	0.01	Stationary
Phillips-Perron Test	-240.783	0.01	Stationary
KPSS Test	0.077	0.1	Stationary
ERS Test	-7.521	N/A	Stationary

check_normality_residuals(ar_011101)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	1.064	0.008	Non-Normal
Jarque-Bera Test	20.178	0	Non-Normal
Lilliefors Test	0.052	0.115	Normal
Skewness	-0.050	N/A	Approximately Symmetric
Kurtosis	4.420	N/A	Leptokurtic (Heavy Tails)

SARIMA Attempt 2 (d=D=1)

In this attempt, we will only focus on SARIMA models with \(d= D=0\).

ARIMA 011010

From the suggestions during differencing, we first fit an ARIMA(0,1,1)(0,1,0)[12]. Based on the ACF and PACF, a SMA 1 (\(Q=1\))should be recalled.

# ARIMA 011010
ar_011010 <- arima(Arrivals_ts, order=c(0,1,1), seasonal=list(order=c(0,1,0), period=12))
ar_011010

Call:
arima(x = Arrivals_ts, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 0), 
    period = 12))

Coefficients:
          ma1
      -0.4833
s.e.   0.0580

sigma^2 estimated as 8.707e+10:  log likelihood = -3181.29,  aic = 6366.59

checkresiduals(ar_011010)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,1)(0,1,0)[12]
Q* = 79.472, df = 23, p-value = 3.865e-08

Model df: 1.   Total lags used: 24

two_plots(ar_011010$residuals)

check_stationarity(ar_011010$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-5.740	0.01	Stationary
Phillips-Perron Test	-258.543	0.01	Stationary
KPSS Test	0.023	0.1	Stationary
ERS Test	-7.880	N/A	Stationary

check_normality_residuals(ar_011010)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	2.302	0	Non-Normal
Jarque-Bera Test	66.311	0	Non-Normal
Lilliefors Test	0.072	0.004	Non-Normal
Skewness	0.030	N/A	Approximately Symmetric
Kurtosis	5.570	N/A	Leptokurtic (Heavy Tails)

ARIMA 011011

Here, we fit an ARIMA(0,1,1)(0,1,1)[12]. Based on the ACF, PACF, and other statistical tests, we have found a model with independent and stationary residuals. However, the residuals are not normal; hence, we cannot trust the forecast intervals.

# ARIMA 011011
ar_011011 <- arima(Arrivals_ts, order=c(0,1,1), 
                   seasonal=list(order=c(0,1,1), period=12))
ar_011011

Call:
arima(x = Arrivals_ts, order = c(0, 1, 1), seasonal = list(order = c(0, 1, 1), 
    period = 12))

Coefficients:
          ma1     sma1
      -0.4600  -0.6602
s.e.   0.0609   0.0510

sigma^2 estimated as 5.947e+10:  log likelihood = -3141.44,  aic = 6288.87

checkresiduals(ar_011011)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,1)(0,1,1)[12]
Q* = 19.923, df = 22, p-value = 0.5879

Model df: 2.   Total lags used: 24

two_plots(ar_011011$residuals)

check_stationarity(ar_011011$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-6.436	0.01	Stationary
Phillips-Perron Test	-243.614	0.01	Stationary
KPSS Test	0.077	0.1	Stationary
ERS Test	-7.635	N/A	Stationary

check_normality_residuals(ar_011011)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	1.374	0.001	Non-Normal
Jarque-Bera Test	19.836	0	Non-Normal
Lilliefors Test	0.073	0.003	Non-Normal
Skewness	-0.010	N/A	Approximately Symmetric
Kurtosis	4.410	N/A	Leptokurtic (Heavy Tails)

iii. Model Validation & Prediction

Model Comparison

Based on the summary table below, the best model appears to be ARIMA(0,1,1)(0,1,1)[12], with the lowest AIC and BIC values. Still, we recall that the residuals for this model are not normal, and thus we should not trust the confidence intervals with much reliability. Yet, the residuals are independent and stationary, which suggests that we should still proceed with forecasting.

ar_sarima_models <- list("auto.arima()"= ar_initial, 
                         "SARIMA(0,1,0)(1,0,0)[12] without drift" = ar_010100,
                         "SARIMA(0,1,1)(1,0,0)[12] without drift" = ar_011100,
                         "SARIMA(0,1,1)(1,0,1)[12] without drift" = ar_011101,
                         "SARIMA(0,1,1)(0,1,0)[12] without drift" = ar_011010, 
                         "SARIMA(0,1,1)(0,1,1)[12] without drift" = ar_011011)
sarima_summary_table(ar_sarima_models)

Model Summary
Model	AIC	BIC
auto.arima()	6320.56	6351.42
SARIMA(0,1,0)(1,0,0)[12] without drift	6764.97	6771.92
SARIMA(0,1,1)(1,0,0)[12] without drift	6716.69	6727.11
SARIMA(0,1,1)(1,0,1)[12] without drift	6647.36	6661.27
SARIMA(0,1,1)(0,1,0)[12] without drift	6366.59	6373.44
SARIMA(0,1,1)(0,1,1)[12] without drift	6288.87	6299.15

check_model(ar_011011)

Cross Validation

After verifying that ARIMA(0,1,1)(0,1,1)[12] is also the best fit for the time period between 2000 and 2018, we used this model to predict the 2019 value and tested against the actual data. The result is astonishing with the fitted values following closely actual values. Moreover, an R^2 value of 0.9707 suggests that 97.07% of the variability in actual data can be explained by our model. We now want to leverage this model to predict the 2020 and 2021 hypothetical scenarios of no COVID.

#Setting up
train_data <- window(Arrivals_ts, end = c(2018, 12))
test_data <- window(Arrivals_ts, start = c(2019, 1), end = c(2019, 12))

forecast_2019 <- sarima.for(train_data, n.ahead = 12, 0, 1, 1, 0, 1, 1, 12, plot =FALSE)

predicted_2019 <- forecast_2019$pred
lower_ci_2019 <- forecast_2019$pred - 1.96 * forecast_2019$se
upper_ci_2019 <- forecast_2019$pred + 1.96 * forecast_2019$se

time_full_2019 <- data.frame(
  Date = c(time(train_data), time(test_data)),
  Actual = c(as.numeric(train_data), as.numeric(test_data)),
  Predicted = c(rep(NA, length(train_data)), predicted_2019),
  Lower_CI = c(rep(NA, length(train_data)), lower_ci_2019),
  Upper_CI = c(rep(NA, length(train_data)), upper_ci_2019))

time_only_2019 <- data.frame(
  Date = as.Date(time(test_data)),
  Actual = test_data,
  Predicted = predicted_2019,
  Lower_CI = lower_ci_2019,
  Upper_CI = upper_ci_2019)

# Plotting only 2019 
ar_pred_2019_sarima <- ggplot(time_only_2019, aes(x = Date)) +
  geom_line(aes(y = Actual / 1000000, color = "Actual")) +
  geom_line(aes(y = Predicted / 1000000, color = "Predicted")) +
  geom_ribbon(aes(ymin = Lower_CI / 1000000, ymax = Upper_CI / 1000000), fill = "grey30", alpha = 0.4) +
  labs(title = "Actual and Predicted Arrivals (2019) with 95% CI",
       y = "Tourist Arrivals (millions)") +
  scale_x_date(date_labels = "%b %Y", date_breaks = "2 months") +
  theme_minimal() +
  scale_color_manual(values = c("Actual" = "darkblue", "Predicted" = "red"))

# Plotting full 
ar_pred_full_sarima <- ggplot(time_full_2019, aes(x = Date)) +
  geom_line(aes(y = Actual / 1000000, color = "Actual")) +
  geom_line(aes(y = Predicted / 1000000, color = "Predicted")) +
  geom_ribbon(aes(ymin = Lower_CI / 1000000, ymax = Upper_CI / 1000000), fill = "grey30", alpha = 0.4) +
  labs(title = "Actual and Predicted Arrivals (2019) with 95% CI",
       y = "Tourist Arrivals (millions)") +
  theme_minimal() +
  scale_color_manual(values = c("Actual" = "darkblue", "Predicted" = "red"))

grid.arrange(ar_pred_full_sarima, ar_pred_2019_sarima, nrow = 2)

# R^Squared 
calculate_r2(as.numeric(predicted_2019), as.numeric(test_data))

R-squared of the model is: 0.9707 

Prediction

Here, we can see that the model predicts monthly arrivals to be relatively stable, exhibiting historical seasonal trends with a slightly upward motion. This makes sense as without COVID-19, the US should have remained its attractiveness as an ideal place for travel, business, and education in 2020 and 2021.

# Run SARIMA forecast for 24 months ahead (2020 and 2021)
ar_sarima_for <- sarima.for(Arrivals_ts, n.ahead = 24, 0,1,1,0,1,1, S = 12, plot = FALSE)

forecasted_values <- ar_sarima_for$pred
lower_bound <- ar_sarima_for$pred - 1.96 * ar_sarima_for$se
upper_bound <- ar_sarima_for$pred + 1.96 * ar_sarima_for$se

time_forecast <- ts(rep(NA,24), start = c(2020, 1), end = c(2021, 12), frequency = 12)

forecast_ts <- ts(forecasted_values, start = c(2020, 1), frequency = 12)
lower_bound_ts <- ts(lower_bound, start = c(2020, 1), frequency = 12)
upper_bound_ts <- ts(upper_bound, start = c(2020, 1), frequency = 12)

combined_ts <- ts(c(Arrivals_ts, forecast_ts), start = start(Arrivals_ts), frequency = 12)

time_full_forecast <- data.frame(
  Date = c(time(Arrivals_ts), time(time_forecast)),
  Actual = c(as.numeric(Arrivals_ts), rep(NA, length(forecast_ts))),
  Predicted = c(rep(NA, length(Arrivals_ts)), forecasted_values),
  Lower_CI = c(rep(NA, length(Arrivals_ts)), lower_bound),
  Upper_CI = c(rep(NA, length(Arrivals_ts)), upper_bound))

time_forecast_2020_2021 <- data.frame(
  Date = as.Date(time(time_forecast)),
  Predicted = forecasted_values,
  Lower_CI = lower_bound,
  Upper_CI = upper_bound)

# Plotting full dataset with forecasts
ar_pred_full_forecast <- ggplot(time_full_forecast, aes(x = Date)) +
  geom_line(aes(y = Actual / 1000000, color = "Actual")) +
  geom_line(aes(y = Predicted / 1000000, color = "Predicted")) +
  geom_ribbon(aes(ymin = Lower_CI / 1000000, ymax = Upper_CI / 1000000), fill = "grey30", alpha = 0.4) +
  labs(title = "Full Data with Forecasted Arrivals (2020-2021) with 95% CI",
       y = "Tourist Arrivals (millions)") +
  theme_minimal() +
  scale_color_manual(values = c("Actual" = "darkblue", "Predicted" = "red"))

# Plotting only 2020 and 2021
ar_pred_2020_2021 <- ggplot(time_forecast_2020_2021, aes(x = Date)) +
  geom_line(aes(y = Predicted / 1000000, color = "Predicted")) +
  geom_ribbon(aes(ymin = Lower_CI / 1000000, ymax = Upper_CI / 1000000), fill = "grey30", alpha = 0.4) +
  labs(title = "Forecasted Tourist Arrivals (2020-2021) with 95% CI",
       y = "Tourist Arrivals (millions)") +
  scale_x_date(date_labels = "%b %Y", date_breaks = "6 months") +
  theme_minimal() +
  scale_color_manual(values = c("Predicted" = "red"))

grid.arrange(ar_pred_full_forecast, ar_pred_2020_2021, nrow = 2)

b. Tourist Spending

i. Achieving Stationarity

Time Series

This time series is definitely not stationary, based on the time series plot, the tests, and the ACF/PACF plots.

diff_plots(Exports_ts/1000, "Tourist Exports, billions USD (2000-2019)")

check_stationarity(Exports_ts)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-3.005	0.15	Non-Stationary
Phillips-Perron Test	-8.956	0.61	Non-Stationary
KPSS Test	4.746	0.01	Non-Stationary
ERS Test	1.206	N/A	Non-Stationary

Auto ARIMA

Auto.arima() suggests an ARIMA(2,1,2)(2,0,0)[12] with drift. From the outset, this model seems to fit the model well, with stationary and independent residuals. Yet, the residuals are not normal, suggesting poor reliability of forecasts. This step is only for reference, and we will attempt to build our only models.

# Auto.arima() - ARIMA(2,1,2)(2,0,0)[12] with drift
ex_initial <- auto.arima(Exports_ts)
ex_initial

Series: Exports_ts 
ARIMA(2,1,2)(2,0,0)[12] with drift 

Coefficients:
          ar1      ar2     ma1     ma2     sar1     sar2    drift
      -0.4818  -0.8840  0.5054  0.9604  -0.2794  -0.1773  47.1808
s.e.   0.0631   0.0615  0.0379  0.0378   0.0685   0.0769  12.5827

sigma^2 = 72552:  log likelihood = -1674.01
AIC=3364.03   AICc=3364.65   BIC=3391.84

ts.plot(cbind(Exports_ts/1000,ex_initial$fitted/1000), 
        col=c("darkblue", "darkred"), ylab = "Exports (billions $)", 
        main = "Original Data with Fitted Values - Auto ARIMA", lwd =1.5)
legend("topleft", legend=c("Original Data", "Fitted Values"), 
       col=c("darkblue", "darkred"), lty=1:1,lwd=2)

checkresiduals(ex_initial)

    Ljung-Box test

data:  Residuals from ARIMA(2,1,2)(2,0,0)[12] with drift
Q* = 17.561, df = 18, p-value = 0.4849

Model df: 6.   Total lags used: 24

two_plots(ex_initial$residuals)

check_stationarity(ex_initial$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-5.506	0.01	Stationary
Phillips-Perron Test	-226.347	0.01	Stationary
KPSS Test	0.297	0.1	Stationary
ERS Test	-6.948	N/A	Stationary

check_normality_residuals(ex_initial)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	1.561	0.001	Non-Normal
Jarque-Bera Test	106.788	0	Non-Normal
Lilliefors Test	0.063	0.021	Non-Normal
Skewness	-0.450	N/A	Approximately Symmetric
Kurtosis	6.140	N/A	Leptokurtic (Heavy Tails)

Differencing

Recall the slow decay in the ACF plot, normal differencing is the first thing we should try. Luckily, it appears to be all that we need to achieve stationarity. Noting the non-zero mean of residuals, we will include a drift term in our models. Based on the ACF and PACF, the best next move appears to be a SMA 1 (\(Q=1\))

ex_diff <- diff(Exports_ts)

diff_plots(ex_diff/1000,"Normally Differenced Time Series")

check_stationarity(ex_diff/1000)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-5.853	0.01	Stationary
Phillips-Perron Test	-233.100	0.01	Stationary
KPSS Test	0.157	0.1	Stationary
ERS Test	-2.648	N/A	Stationary

ii. Model Fitting

SARIMA Attempt 1

From the suggestions during differencing, we first fit an ARIMA(0,1,0)(0,0,1)[12]. Luckily, we immediately arrived at a good model with residuals showing stationarity and independence. However, the residuals are not normal; hence, we cannot trust the forecast intervals.

ex_010001 <- Arima(Exports_ts, order=c(0,1,0), seasonal=list(order=c(0,0,1), period=12), include.drift=TRUE)
ex_010001

Series: Exports_ts 
ARIMA(0,1,0)(0,0,1)[12] with drift 

Coefficients:
         sma1    drift
      -0.2678  46.2844
s.e.   0.0725  13.1017

sigma^2 = 74133:  log likelihood = -1678.6
AIC=3363.19   AICc=3363.29   BIC=3373.62

checkresiduals(ex_010001)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,0)(0,0,1)[12] with drift
Q* = 21.046, df = 23, p-value = 0.5783

Model df: 1.   Total lags used: 24

two_plots(ex_010001$residuals)

check_stationarity(ex_010001$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-5.669	0.01	Stationary
Phillips-Perron Test	-227.123	0.01	Stationary
KPSS Test	0.267	0.1	Stationary
ERS Test	-6.824	N/A	Stationary

check_normality_residuals(ex_010001)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	1.577	0	Non-Normal
Jarque-Bera Test	115.892	0	Non-Normal
Lilliefors Test	0.083	0	Non-Normal
Skewness	-0.430	N/A	Approximately Symmetric
Kurtosis	6.290	N/A	Leptokurtic (Heavy Tails)

SARIMA Attempt 2

ARIMA 010200 w/ drift

From the ACF anf PACF plots during differencing, someone may suggest an SAR2 as well, i.e. (\(P=2\)). Hence, we fit an ARIMA(0,1,0)(2,0,0)[12]. Here, we also achieved stationarity and independence for residuals, though not normality. With two solid models, we wanted to see what would happen if we include more terms.

ex_010200 <- Arima(Exports_ts, order=c(0,1,0), seasonal=list(order=c(2,0,0), period=12), include.drift=TRUE)
ex_010200

Series: Exports_ts 
ARIMA(0,1,0)(2,0,0)[12] with drift 

Coefficients:
         sar1     sar2    drift
      -0.2526  -0.1882  47.2346
s.e.   0.0665   0.0754  12.3734

sigma^2 = 73283:  log likelihood = -1676.98
AIC=3361.95   AICc=3362.12   BIC=3375.86

checkresiduals(ex_010200)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,0)(2,0,0)[12] with drift
Q* = 22.37, df = 22, p-value = 0.438

Model df: 2.   Total lags used: 24

two_plots(ex_010200$residuals)

check_stationarity(ex_010200$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-5.651	0.01	Stationary
Phillips-Perron Test	-222.324	0.01	Stationary
KPSS Test	0.295	0.1	Stationary
ERS Test	-6.752	N/A	Stationary

check_normality_residuals(ex_010200)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	1.557	0.001	Non-Normal
Jarque-Bera Test	103.015	0	Non-Normal
Lilliefors Test	0.077	0.001	Non-Normal
Skewness	-0.410	N/A	Approximately Symmetric
Kurtosis	6.100	N/A	Leptokurtic (Heavy Tails)

ARIMA 011200 w/ drift

ex_011200 <- Arima(Exports_ts, order=c(0,1,1), seasonal=list(order=c(2,0,0), period=12), include.drift=TRUE)
ex_011200

Series: Exports_ts 
ARIMA(0,1,1)(2,0,0)[12] with drift 

Coefficients:
         ma1     sar1     sar2    drift
      0.0786  -0.2635  -0.1988  47.4203
s.e.  0.0602   0.0669   0.0752  13.1025

sigma^2 = 73033:  log likelihood = -1676.13
AIC=3362.27   AICc=3362.52   BIC=3379.65

checkresiduals(ex_011200)

    Ljung-Box test

data:  Residuals from ARIMA(0,1,1)(2,0,0)[12] with drift
Q* = 20.941, df = 21, p-value = 0.4626

Model df: 3.   Total lags used: 24

two_plots(ex_011200$residuals)

check_stationarity(ex_011200$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-5.667	0.01	Stationary
Phillips-Perron Test	-243.946	0.01	Stationary
KPSS Test	0.294	0.1	Stationary
ERS Test	-6.871	N/A	Stationary

check_normality_residuals(ex_011200)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	1.527	0.001	Non-Normal
Jarque-Bera Test	84.386	0	Non-Normal
Lilliefors Test	0.070	0.007	Non-Normal
Skewness	-0.320	N/A	Approximately Symmetric
Kurtosis	5.830	N/A	Leptokurtic (Heavy Tails)

ARIMA 111200 w/ drift

ex_111200 <- Arima(Exports_ts, order=c(1,1,1), seasonal=list(order=c(2,0,0), period=12), include.drift=TRUE)
ex_111200

Series: Exports_ts 
ARIMA(1,1,1)(2,0,0)[12] with drift 

Coefficients:
         ar1      ma1     sar1     sar2    drift
      0.4505  -0.3556  -0.2654  -0.1967  47.5588
s.e.  0.3912   0.4082   0.0669   0.0749  14.2048

sigma^2 = 73055:  log likelihood = -1675.66
AIC=3363.32   AICc=3363.68   BIC=3384.18

checkresiduals(ex_111200)

    Ljung-Box test

data:  Residuals from ARIMA(1,1,1)(2,0,0)[12] with drift
Q* = 20.388, df = 20, p-value = 0.4339

Model df: 4.   Total lags used: 24

two_plots(ex_111200$residuals)

check_stationarity(ex_111200$residuals)

Stationarity Test Summary
Test	Test.Statistic	P.Value	Interpretation
ADF Test	-5.863	0.01	Stationary
Phillips-Perron Test	-240.199	0.01	Stationary
KPSS Test	0.272	0.1	Stationary
ERS Test	-7.128	N/A	Stationary

check_normality_residuals(ex_111200)

Normality Test Summary
Test	Test.Statistic	P.Value	Interpretation
Anderson-Darling Test	1.663	0	Non-Normal
Jarque-Bera Test	91.445	0	Non-Normal
Lilliefors Test	0.066	0.013	Non-Normal
Skewness	-0.290	N/A	Approximately Symmetric
Kurtosis	5.970	N/A	Leptokurtic (Heavy Tails)

iii. Model Validation & Prediction

Model Comparison

Of all five models, our first two models appear to be the best fits to the data. However, while the second one (SARIMA(0,1,0)(2,0,0)[12] with drift) is more favorable in terms of AIC, the first one (SARIMA(0,1,0)(0,0,1)[12] with drift) is superior in terms of BIC. Now, as we know that the time series is seasonally adjusted, which means that seasonality should not play too much of a role in our prediction; hence, the notion that tourist exports from two years ago may still have an impact on that of this current month seems unplausible. Therefore, we believe that the first model, SARIMA(0,1,0)(0,0,1)[12] with drift, should be deemed as best fit. Again, as the residuals are not normal, we should not blindly trust our forecast intervals.

ex_sarima_models <- list("auto.arima() "= ex_initial,
                         "SARIMA(0,1,0)(0,0,1)[12] with drift" = ex_010001,
                         "SARIMA(0,1,0)(2,0,0)[12] with drift" = ex_010200,
                         "SARIMA(0,1,1)(2,0,0)[12] with drift" = ex_011200,
                         "SARIMA(1,1,1)(2,0,0)[12] with drift" = ex_111200)
sarima_summary_table(ex_sarima_models)

Model Summary
Model	AIC	BIC
auto.arima()	3364.03	3391.84
SARIMA(0,1,0)(0,0,1)[12] with drift	3363.19	3373.62
SARIMA(0,1,0)(2,0,0)[12] with drift	3361.95	3375.86
SARIMA(0,1,1)(2,0,0)[12] with drift	3362.27	3379.65
SARIMA(1,1,1)(2,0,0)[12] with drift	3363.32	3384.18

check_model(ex_010001)

Cross Validation

After verifying that ARIMA(0,1,0)(0,0,1)[12] with drift is also the best fit for the time period between 2000 and 2018, we used this model to predict the 2019 value and tested against the actual data. While the fitted values do not follow the actual data too closely, the actual data falls under the 95% interval, which is a good sign that our model has some predictive power. However, a negative R^2 value of -1.6157 is concernning, and needs to be taken seriously later. We now want to leverage this model to predict the 2020 and 2021 hypothetical scenarios of no COVID.

# Setting up
train_data_exports <- window(Exports_ts, end = c(2018, 12))
test_data_exports <- window(Exports_ts, start = c(2019, 1), end = c(2019, 12))

forecast_exports <- sarima.for(train_data_exports, n.ahead = 12, 0, 1, 0, 0, 0, 1, S = 12, plot = FALSE, no.constant = FALSE)

predicted_exports <- forecast_exports$pred
lower_ci_exports <- forecast_exports$pred - 1.96 * forecast_exports$se
upper_ci_exports <- forecast_exports$pred + 1.96 * forecast_exports$se

time_full_exports <- data.frame(
  Date = c(time(train_data_exports), time(test_data_exports)),
  Actual = c(as.numeric(train_data_exports), as.numeric(test_data_exports)),
  Predicted = c(rep(NA, length(train_data_exports)), predicted_exports),
  Lower_CI = c(rep(NA, length(train_data_exports)), lower_ci_exports),
  Upper_CI = c(rep(NA, length(train_data_exports)), upper_ci_exports))

time_only_exports <- data.frame(
  Date = as.Date(time(test_data_exports)),
  Actual = as.numeric(test_data_exports),
  Predicted = predicted_exports,
  Lower_CI = lower_ci_exports,
  Upper_CI = upper_ci_exports)

# Plotting only 2020-2021 forecast
ar_pred_exports_2020_2021 <- ggplot(time_only_exports, aes(x = Date)) +
  geom_line(aes(y = Actual / 1000, color = "Actual")) +
  geom_line(aes(y = Predicted / 1000, color = "Predicted")) +
  geom_ribbon(aes(ymin = Lower_CI / 1000, ymax = Upper_CI / 1000), fill = "grey30", alpha = 0.4) +
  labs(title = "Actual and Predicted Exports (2020-2021) with 95% CI",
       y = "Exports (billions $)") +
  scale_x_date(date_labels = "%b %Y", date_breaks = "2 months") +
  theme_minimal() +
  scale_color_manual(values = c("Actual" = "darkblue", "Predicted" = "red"))

# Plotting full dataset with forecast
ar_pred_full_exports <- ggplot(time_full_exports, aes(x = Date)) +
  geom_line(aes(y = Actual / 1000, color = "Actual")) +
  geom_line(aes(y = Predicted / 1000, color = "Predicted")) +
  geom_ribbon(aes(ymin = Lower_CI / 1000, ymax = Upper_CI / 1000), fill = "grey30", alpha = 0.4) +
  labs(title = "Full Data with Predicted Exports (2020-2021) and 95% CI",
       y = "Exports (billions $)") +
  theme_minimal() +
  scale_color_manual(values = c("Actual" = "darkblue", "Predicted" = "red"))

grid.arrange(ar_pred_full_exports, ar_pred_exports_2020_2021, nrow = 2)

calculate_r2(as.numeric(predicted_exports), as.numeric(test_data_exports))

R-squared of the model is: -1.6157 

Prediction

Here, we can see that the model predicts monthly exports to be relatively stable, with a slightly upward motion. This makes sense as without COVID-19, the US should have remained its attractiveness as an ideal place for travel, business, and education, and continued to grow their tourism exports in 2020 and 2021.

# Run SARIMA forecast
ex_sarima_for <- sarima.for(Exports_ts, n.ahead = 24, 0, 1, 0, 0, 0, 1, S = 12, 
                            plot = FALSE, no.constant = FALSE)

forecasted_values <- ex_sarima_for$pred
lower_bound <- ex_sarima_for$pred - 1.96 * ex_sarima_for$se
upper_bound <- ex_sarima_for$pred + 1.96 * ex_sarima_for$se

time_forecast <- ts(rep(NA, 24), start = c(2020, 1), end = c(2021, 12), frequency = 12)
forecast_ts <- ts(forecasted_values, start = c(2020, 1), frequency = 12)
lower_bound_ts <- ts(lower_bound, start = c(2020, 1), frequency = 12)
upper_bound_ts <- ts(upper_bound, start = c(2020, 1), frequency = 12)

combined_ts <- ts(c(Exports_ts, forecast_ts), start = start(Exports_ts), frequency = 12)

time_full_forecast <- data.frame(
  Date = c(time(Exports_ts), time(time_forecast)),
  Actual = c(as.numeric(Exports_ts), rep(NA, length(forecast_ts))),
  Predicted = c(rep(NA, length(Exports_ts)), forecasted_values),
  Lower_CI = c(rep(NA, length(Exports_ts)), lower_bound),
  Upper_CI = c(rep(NA, length(Exports_ts)), upper_bound))


time_forecast_2020_2021 <- data.frame(
  Date = as.Date(time(time_forecast)),
  Predicted = forecasted_values,
  Lower_CI = lower_bound,
  Upper_CI = upper_bound)

# Plotting full dataset with forecasts
ex_pred_full_forecast <- ggplot(time_full_forecast, aes(x = Date)) +
  geom_line(aes(y = Actual / 1000, color = "Actual")) +
  geom_line(aes(y = Predicted / 1000, color = "Predicted")) +
  geom_ribbon(aes(ymin = Lower_CI / 1000, ymax = Upper_CI / 1000), fill = "grey30", alpha = 0.4) +
  labs(title = "Full Data with Forecasted Exports (2020-2021) with 95% CI",
       y = "Exports (billions $)") +
  theme_minimal() +
  scale_color_manual(values = c("Actual" = "darkblue", "Predicted" = "red"))

# Plotting only 2020–2021 forecast
ex_pred_2020_2021 <- ggplot(time_forecast_2020_2021, aes(x = Date)) +
  geom_line(aes(y = Predicted / 1000, color = "Predicted")) +
  geom_ribbon(aes(ymin = Lower_CI / 1000, ymax = Upper_CI / 1000), fill = "grey30", alpha = 0.4) +
  labs(title = "Forecasted Exports (2020-2021) with 95% CI",
       y = "Exports (billions $)") +
  scale_x_date(date_labels = "%b %Y", date_breaks = "3 months") +
  theme_minimal() +
  scale_color_manual(values = c("Predicted" = "red"))

grid.arrange(ex_pred_full_forecast, ex_pred_2020_2021, nrow = 2)

4. Cross Correlation

Based on the two CCF plots, it can be seen that the changes in this month’s tourist arrivals can be attributed to changes in tourist exports 10 and 11 months earlier. As these relationships seem to be odd, future analysis will dive deeper into these analyses.

exports_model <- sarima(Exports_ts, p=0, d=1, q=0, P=0, D=0, Q=1, S=12,details=FALSE) 

<><><><><><><><><><><><><><>
 
Coefficients: 
         Estimate      SE t.value p.value
sma1      -0.2678  0.0725 -3.6943   3e-04
constant  46.2844 13.1017  3.5327   5e-04

sigma^2 estimated as 73512.38 on 237 degrees of freedom 
 
AIC = 14.07193  AICc = 14.07214  BIC = 14.11557 
 

Exports_w <- resid(exports_model$fit) 
phi <- as.vector(exports_model$fit$coef) 

Arrivals_f1 <- stats::filter(diff(Arrivals_sa_saar), method ="convolution", filter = c(rep(0,11),-phi[1]), sides = 1) 

both1 <- ts.intersect(Exports_w, Arrivals_f1) 
Exports_w_aligned1 <- both1[,1] 
Arrivals_f_aligned1 <- both1[,2] 


Arrivals_f2 <- stats::filter(diff(Arrivals_sa_stl), method ="convolution", filter = c(rep(0,11),-phi[1]), sides = 1) 

both2 <- ts.intersect(Exports_w, Arrivals_f2) 
Exports_w_aligned2 <- both2[,1] 
Arrivals_f_aligned2 <- both2[,2] 



ccf(Exports_w_aligned1, Arrivals_f_aligned1, main = "CCF of Whitened Exports and Filtered Arrivals - SAAR", na.action = na.pass)

ccf(Exports_w_aligned2, Arrivals_f_aligned2, main = "CCF of Whitened Exports and Filtered Arrivals - STL", na.action = na.pass)

IV. Closing

This project has been a meaningful learning experience that allowed me to consolidate and expand my understanding of time series analysis. By working independently through each phase, I gained a deeper insight into structuring a comprehensive approach to building predictive models from scratch, applying a series of tests and functions to evaluate model suitability. Additionally, I developed custom functions to streamline this analysis, making it easier to interpret results and refine the model-building-and-validation process. For future research, I would be interested in exploring GARCH models to address potential non-stationary variance in residuals and investigating possible relationships between tourist arrivals and exports, a connection I was unable to examine fully due to time constraints. Overall, this project has given me a clearer sense of the concepts and techniques essential to building robust time series models and has shown me areas to revisit in my coursework to further strengthen my skills.

V. Appendix & References

1. Statistical Tests & Methods

Below is a brief introduction to each statistical test and method utilized in this analysis. Sources of information that helped in compiling these summaries are listed in the References section.

ADF (Augmented Dickey-Fuller) Test

The ADF test checks if a time series is stationary by testing for a unit root. The null hypothesis (H_0) is that the series has a unit root (non-stationary), while the alternative hypothesis (H_1) is that the series is stationary. A significant p-value (typically < 0.05) leads to rejecting H_0, indicating that the series is likely stationary.

Phillips-Perron (PP) Test

Like the ADF test, the PP test evaluates whether a time series is stationary, accounting for potential serial correlation and heteroskedasticity in errors. Here, H_0 states that the series has a unit root (non-stationary), while rejecting H_0 indicates stationarity. A significant result supports stationarity.

KPSS (Kwiatkowski-Phillips-Schmidt-Shin) Test

The KPSS test approaches stationarity testing differently, assuming the series is stationary under H_0. If the p-value is small, H_0 is rejected, suggesting the series is non-stationary. This complements ADF and PP tests by testing stationarity as the null.

ERS (Elliott-Rothenberg-Stock) Test

The ERS test also assesses unit roots, particularly when seeking stationarity in the presence of trend and constant terms. H_0 states the series has a unit root (non-stationary), and rejecting H_0 suggests stationarity. This test is often more sensitive than the ADF test.

ACF (Autocorrelation Function)

The ACF measures the correlation of a series with its own lagged values across different time lags. Significant spikes in the ACF plot at specific lags can indicate autocorrelation, providing insight into patterns and cycles within the data, helping identify moving average terms.

PACF (Partial Autocorrelation Function)

PACF measures the correlation between a series and its lagged values after removing correlations attributed to intermediate lags. Significant spikes in the PACF help determine the order of autoregressive terms in models, particularly in ARIMA modeling.

Adjusted R^2

Adjusted R^2 measures the proportion of variance explained by the model, adjusting for the number of predictors to avoid overfitting. A higher Adjusted R^2 indicates a better-fitting model without adding unnecessary complexity.

AIC (Akaike Information Criterion)

AIC is used to compare models, considering both fit and complexity. A lower AIC suggests a better model, as it balances model fit and the number of parameters. It penalizes excessive complexity to favor more parsimonious models.

BIC (Bayesian Information Criterion)

Similar to AIC but with a larger penalty for additional parameters, BIC also aids in model selection. Lower BIC values indicate a preferable model. BIC is more conservative than AIC, making it useful when simplicity is prioritized.

Anderson-Darling (AD) Test

The AD test is a normality test that gives more weight to the tails of the distribution. Under H_0, the data is normally distributed; rejection indicates significant deviations from normality, especially in the distribution tails.

Jarque-Bera (JB) Test

The JB test evaluates whether skewness and kurtosis of a dataset match a normal distribution. The null hypothesis H_0 is that the data is normally distributed; rejection suggests the presence of non-normality, often due to skewed or heavy-tailed distributions.

Lilliefors Test

A variation of the Kolmogorov-Smirnov test for normality, used when population parameters are unknown. H_0 assumes normality, and a significant p-value indicates a non-normal distribution, especially useful for unknown mean and variance.

Skewness

Skewness measures the asymmetry of data around its mean. A skewness close to zero suggests symmetry, while positive or negative values indicate right or left skew, respectively. Excessive skew may suggest transformations are necessary for normality.

Kurtosis

Kurtosis quantifies the “tailedness” of data, with values near 3 indicating normality (mesokurtic). Higher kurtosis (leptokurtic) suggests heavy tails, while lower kurtosis (platykurtic) indicates light tails.

2. References

“Augmented Dickey Fuller Test (ADF Test) – Must Read Guide.” 2019. Machine Learning Plus (blog). November 2, 2019. https://www.machinelearningplus.com/time-series/augmented-dickey-fuller-test/.

Bevans, Rebecca. 2020. “Akaike Information Criterion | When & How to Use It (Example).” Scribbr. March 26, 2020. https://www.scribbr.com/statistics/akaike-information-criterion.

Datalab, Analyttica. 2019. “What Is Bayesian Information Criterion (BIC)?” Medium. January 16, 2019. https://medium.com/@analyttica/what-is-bayesian-information-criterion-bic-b3396a894be6.

Frost, Jim. 2021. “Autocorrelation and Partial Autocorrelation in Time Series Data.” Statistics by Jim. May 17, 2021. https://statisticsbyjim.com/time-series/autocorrelation-partial-autocorrelation/.

Glen, Stephanie. 2016. “Jarque-Bera Test.” Statistics How To. May 8, 2016. https://www.statisticshowto.com/jarque-bera-test/.

Gonzalez-Barrera, Ana. 2020. “After Surging in 2019, Migrant Apprehensions at U.S.-Mexico Border Fell Sharply in Fiscal 2020.” Pew Research Center. November 4, 2020. https://www.pewresearch.org/short-reads/2020/11/04/after-surging-in-2019-migrant-apprehensions-at-u-s-mexico-border-fell-sharply-in-fiscal-2020-2/.

Hinton, Thomas. 2024. “Infographic: Travel and Tourism Drive close to 10% of the US Economy.” Statista Daily Data. Statista. June 20, 2024. https://www.statista.com/chart/32465/travel-and-tourism-contribution-to-gdp-in-the-us/.

“I-94 Arrivals Program.” 2024. International Trade Administration | Trade.gov. 2024. https://www.trade.gov/i-94-arrivals-program.

IBM. 2023. “Adjusted R Squared.” IBM. January 3, 2023. https://www.ibm.com/docs/en/cognos-analytics/11.1.0?topic=terms-adjusted-r-squared.

“International Travel Receipts and Payments Program.” 2024. International Trade Administration | Trade.gov. 2024. https://www.trade.gov/international-travel-receipts-and-payments-program.

Levin, Adam G. 2023. “U.S. Tourism: Economic Impacts and Pandemic Recovery.” Congressional Research Service. https://crsreports.congress.gov/product/pdf/R/R47857.

Malato, Gianluca. 2023. “An Introduction to the Shapiro-Wilk Test for Normality | Built In.” Builtin. 2023. https://builtin.com/data-science/shapiro-wilk-test.

Olivo, Natalie. 2024. “Medium.” Medium. 2024. https://codeburst.io/cross-validation-calculating-r.

“Pp.test Function - RDocumentation.” n.d. R-Project.org. https://www.rdocumentation.org/packages/aTSA/versions/3.1.2.1/topics/pp.test.

“R: Elliott, Rothenberg and Stock Unit Root Test.” 2024. R-Project.org. 2024. https://search.r-project.org/CRAN/refmans/urca/html/ur.ers.html.

“R: Kwiatkowski-Phillips-Schmidt-Shin Test.” n.d. RDocumentation. https://search.r-project.org/CRAN/refmans/aTSA/html/kpss.test.html.

Scott. 2020. “Parameters Primer: ACF (Autocorrelation Function) - Michigan Metrology.” Michigan Metrology. August 16, 2020. https://michmet.com/parameters-primer-acf-autocorrelation-function.

“Seasonal Adjustments and Deseasonalising Data.” 2024. Mathspace.co. 2024. https://mathspace.co/textbooks/syllabuses/Syllabus-845/topics/Topic-18558/subtopics/Subtopic-251579/?activeTab=theory.

“Skewness | R Tutorial.” 2020. R-Tutor. 2020. https://www.r-tutor.com/elementary-statistics/numerical-measures/skewness.

Stephanie. 2016. “Lilliefors Test for Normality & Exponential Distributions.” Statistics How To. March 8, 2016. https://www.statisticshowto.com/lilliefors-test/.

“The Anderson-Darling Statistic.” n.d. Minitab. https://support.minitab.com/en-us/minitab/help-and-how-to/statistics/basic-statistics/supporting-topics/normality/the-anderson-darling-statistic/.

Zach. 2020. “How to Calculate Skewness & Kurtosis in R.” Statology. October 23, 2020. https://www.statology.org/skewness-kurtosis-in-r/.

Zaheer, Aima. 2023. “25 States with Highest Tourism Revenue in the US.” Yahoo Finance. November 29, 2023. https://finance.yahoo.com/news/25-states-highest-tourism-revenue-151108366.html.